A risk model based on pyroptosis subtypes predicts tumor immune microenvironment and guides chemotherapy and immunotherapy in bladder cancer

Although immunotherapy has revolutionized bladder cancer (BLCA) therapy, only few patients demonstrate durable clinical benefits due to the heterogeneity. Emerging evidence has linked pyroptosis to shaping tumor microenvironment (TME) and predicting therapy response. However, the relationship between pyroptosis and immunotherapy response in BLCA remains elusive. In this study, we performed a comprehensive bioinformatic analysis to dissect the role of pyroptosis in BLCA. Differentially expressed pyroptosis-related genes (DEPRGs) between tumor and normal tissues were identified using publicly available datasets. Kaplan–Meier analysis was performed to screen for DEPRGs associated with survival. Consensus clustering was used for BLCA subtyping. TME characteristics were evaluated by CIBERSORT, ESTIMATE and immune checkpoint genes (ICGs). Following univariate COX regression and LASSO analyses with pyroptosis-related DEGs, the risk model and nomogram were constructed with TCGA dataset and validated in the GEO dataset. Furthermore, therapeutic responses in high- and low-risk groups were compared using TIDE and GDSC databases. Two pyroptosis-related subtypes (Cluster 1 and 2) were identified based on expression patterns of GSDMA and CHMP4C. Bioinformatic analyses showed that cluster 1 had poor survival, more M0/M1/M2 macrophages, higher immune/stromal/ESTIMATE scores, and higher expression levels of ICGs. A 15-gene signature for predicting prognosis could classify patients into high- and low-risk groups. Furthermore, the correlation of risk scores with TIDE score and IC50 showed that patients in low-risk group were more sensitive to immunotherapy, whereas patients in high-risk group could better benefit from chemotherapy. Our study identified two novel pyroptosis-related subtypes and constructed a risk model, which can predict the prognosis, improve our understanding the role of PRGs in BLCA, and guide chemotherapy and immunotherapy.

BLCA patients to novel therapeutic agents are limited, predisposing them to inferior clinical outcomes 6 . Although intrinsic molecular subtypes (luminal/basal) have been used for prognostication and therapy response prediction, the clinical value of such subtypes in the prognosis of BLCA remains controversial 2 . Recent studies have highlighted the important roles of the molecular characteristics of patients in improving the therapeutic effects for advanced BLCA 7,8 . Therefore, there is an urgent need to develop an effective gene signature-based risk model in order to guide clinical treatments, especially for chemotherapy and immunotherapy.
Pyroptosis, a newly discovered form of programmed cell death, is initially described in myeloid cells infected by pathogens in 1992 9 . In recent years, the biggest breakthrough in this field has been the elucidation of the roles of gasdermin family as executors of pyroptosis 10,11 . Gasdermins can be cleaved by caspase-1/4/5/11 to release their N-terminal domains, which perforate cell membranes to induce membrane rupture and the release of cytosolic contents into the extracellular environment, amplifying the local or systemic inflammatory response and activating the cellular pyroptosis 10,12 . In the past decade, emerging evidence has demonstrated that pyroptosis could also play important roles in the cancer pathophysiology, regulation of tumor immune microenvironment and modulating the efficacy of immunotherapy and chemotherapy 13 . Zhang et al. found that gasdermin E (GSDME) could function as tumor suppressor gene by activating pyroptosis, which was cleaved by caspase 3 and granzyme B, enhancing macrophage-mediated anti-tumor immunity 14 . Wang et al. 15 also used a biorthogonal system to further verify the robust pyroptosis-induced anti-tumor immunity in a mouse model, which could synergize with cancer immunotherapy using checkpoint blockade. In addition, many chemotherapeutic agents could induce caspase-3/GSDME-dependent pyroptosis in normal cells, which might explain the side effects of chemotherapy [16][17][18] . However, a recent study found that pyroptosis occurred in tumor cells located in the central hypoxia region, contributing to chronic tumor necrosis and a dampened anti-tumor immunity, which in turn promoted tumor progression 19 . Thus, the above studies have shown that the role of pyroptosis in anti-tumor or tumor-promoting remains controversial 20,21 . Extensive research is still needed to clarify the role of pyroptosis in the development of cancers, especially in the understudied BLCA.
In this study, we aim to build a pyroptosis-related prognostic model to predict the prognosis, evaluate the tumor infiltrating immune landscape, and inform chemotherapy and immunotherapy response in BLCA. First, 414 BLCA samples from The Cancer Genome Atlas (TCGA) database were divided into two subtypes based on pyroptosis-related genes (PRGs). A total of 953 differentially expressed genes (DEGs) between two subtypes were identified and then used to constructed a prognostic model. Univariate and multivariate COX regression as well as the least absolute shrinkage and selection operation (LASSO) were used to screen the prognosis-related DEGs, and the risk model was constructed and validated in the Gene Expression Omnibus (GEO) cohort. Finally, a nomogram which combined with clinical information and prognostic risk score was established to comprehensively predict BLCA patients' overall survival (OS). Collectively, our study showed the prognostic value of pyrotosis-related risk model, which had the potential to predict clinical outcomes and responses to chemotherapy and immunotherapy for BLCA.

Materials and methods
Data collection and processing. All datasets used in this study were collected from the public available databases, including TCGA and GEO. The RNA-seq expression data and matched clinical information of 414 BLCA patients were downloaded from TCGA database (https:// portal. gdc. cancer. gov/), and were used as the training dataset. GSE13507 dataset, including RNA expression data and clinical information of 246 BLCA patients, was sourced from GEO database (https:// www. ncbi. nlm. nih. gov/ geo/), and was used as external validation dataset. The RNA expression data of GSE13507 dataset, generated by Illumina human-6 v2.0 expression beadchip, was converted to transcripts per kilobase million (TPM) format using R package limma for downstream analyses 22 . 57 pyroptosis-related genes (PRGs) were obtained from REACTOME database (https:// react ome. org/), Gene Ontology (GO) database (http:// geneo ntolo gy. org), and published articles 23 (Table S1). The whole analysis pipeline was shown in the flow chart ( Figure S1).

Identification of survival-related differentially expressed PRGs in BLCA. DESeq2 package in R
software was utilized to screen the differentially expressed genes (DEGs) between 414 BLCA and 19 normal bladder mucosa samples in the training dataset with the criteria of |log 2 (Fold change (FC))|> 1 and adjusted P value < 0.05 24 . Then, the differentially expressed PRGs (DEPRGs) were obtained by interesting DEGs with 57 PRGs. Next, BLCA patients in the training dataset were divided into low-and high-expression groups according to the median expression level of each DEPRG. The survival between low-and high-expression group was analyzed and compared by Kaplan-Meier curves and log-rank test analysis. DEPRGs were considered to be associated with the overall survival of BLCA patients when there was significant survival difference between low-and high-expression groups (P value < 0.05).

Consensus clustering of BLCA patients based on survival-related DEPRGs. Consensus clustering
was performed by ConsensusClusterPlus package to divide the BLCA patients into distinct pyroptosis-related subgroups based on the expression of survival-related DEPRGs 25 . The number of the clusters were determined by k value, and the stability of classification was guaranteed by performing 1,000 times repetitions. t-distributed stochastic neighbor embedding (t-SNE) was further conducted to assess the performance of consensus clustering. Furthermore, to characterize the subgroups, we (1) analyzed and compared the survival of patients in different subgroups by Kaplan-Meier curves and log-rank test, (2) analyzed and compared the distribution of clinical features (age, gender and TNM stage) in different groups, (3) performed the gene set enrichment analysis (GSEA) with R package GSVA to explore the significantly enriched biological processes and pathways between different groups using gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) 26 reference Analysis of immunotherapeutic response and chemotherapeutic sensitivity. Tumor Immune Dysfunction and Exclusion (TIDE) score of each sample in high-and low-risk groups was calculated using TIDE database (http:// tide. dfci. harva rd. edu/) 31 . Difference of TIDE prediction score between high-and low-risk groups was analyzed by Wilcoxon rank sum test. High TIDE score indicated poor outcome of immune checkpoint blockade (ICB) and short survival time after treatment. In addition, based on Genomics of Drug Sensitivity in Cancer (GDSC) online database, the half maximal inhibitory concentration (IC 50 ) value of cisplatin for each patient was calculated 32 , and the differences between high-and low-risk groups were also determined by Wilcoxon rank sum test.
Statistical analysis. Continuous variables were compared between two subgroups using Student's t test for normally distributed variables, and Wilcoxon rank sum test was used to estimate the statistical significance of non-normally distributed variables. Fisher's test or Chi-square test was used for categorical variables. The Kaplan-Meier method with a two-sided log rank test was used to compare the OS of patients between subgroups. Univariate and multivariate Cox proportional hazards regression models were used to assess the independent prognostic value of the risk model. Time-dependent ROC curve analyses were used to estimate the accuracy of the predictive and prognostic value of the risk gene signature model. All statistical analyses were performed with R software (v4.1.0). P value less than 0.05 was considered to be statistically significant.

Results
Pyroptosis-related GSDMA and CHMP4C were related to the development and progression of BLCA. A total of 2,794 DEGs, including 1,509 upregulated and 1,285 downregulated genes, were identified between BLCA and normal tissues (Table S2 and Fig. 1A). By intersecting 2,794 DEGs with 57 PRGs, nine DEPRGs, including IL6, GSDMA, CASP6, ZBP1, CHMP4C, IL1A, AIM2, NLRP2, and NLRP7, were obtained (Fig. 1B). The expression levels of DEPRGs were displayed in the heatmap, indicating that GSDMA, CASP6, ZBP1, CHMP4C, IL1A, AIM2, NLRP2, and NLRP7 were significantly elevated, whereas the expression of IL6 was markedly reduced in BLCA samples compared with those in controls (Fig. 1C). Further Kaplan-Meier analysis showed that BLCA patients with higher expression of GSDMA and CHMP4C had better overall survival (Fig. 2), indicating that pyroptosis-related GSDMA and CHMP4C had close relationship with BLCA survival. To further explore the expression patterns of PRGs in BLCA, a consensus clustering algorithm was used to categorize the BLCA patients based on the expression patterns of GSDMA and CHMP4C. Our results showed that k = 2 was the optimal parameter for classifying the entire TCGA training set into cluster 1 (n = 209) and cluster 2 (n = 199) (Fig. 3A,B). The t-SNE analysis showed the distinct differences of samples between these two clusters (Fig. 3C), further demonstrating the reliability of consensus clustering. Interestingly, we observed a better survival of patients in cluster 2 (log-rank test, P = 0.015, Fig. 3D). Furthermore, the comparisons of clinicopathological features between these two clusters revealed that the distributions of patients at T and M stage were significantly different (Fig. 3E), suggesting that pyroptosis-related GSDMA and CHMP4C had close relationship with the development and progression of BLCA.

Two subtypes had distinct characteristics of tumor microenvironment. Numerous reports have
shown the close relationship between pyroptosis and TME [33][34][35] . Thus, we further investigated the characteristics of TME in the two pyroptosis-related clusters. We first analyzed the immune landscape (22 human immune cells) of the two clusters using the CIBERSORT algorithm (Fig. 4A), and found that the infiltration levels of activated dendritic cells (DCs) were significantly higher in cluster 2 than those in cluster 1, while M0/M1/M2 macrophages had significantly higher infiltration in cluster 1 compared to cluster 2 (Fig. 4B). Moreover, we also calculated the TME score (stromal score, immune score, and ESTIMATE score) of the two clusters using the ESTIMATE algorithm. The results demonstrated that cluster 1 had markedly higher immune, stromal, and ESTIMATE scores than cluster 2 (Fig. 4C). In addition, the expression levels of 45 ICGs were also significantly different between these two clusters, which had much higher expression levels in the majority of ICGs in cluster www.nature.com/scientificreports/ 1 compared with those in cluster 2 (Fig. 4D). These above results indicated that clustering subtypes based on the gene expression profiles of pyroptosis-related GSDMA and CHMP4C were closely associated with TME in BLCA.

Identification of DEGs between pyroptosis-related clusters and functional enrichment analysis.
TO explore the potential biological functions of the two pyroptosis patterns, we identified a total of 953 pyroptosis-related DEGs using the R package "DESeq2" (Table S3 and Fig. 5A) and performed functional enrichment analysis. These DEGs were significantly enriched in 1049 biological processes (BPs), 84 cellular components (CCs), 87 molecular functions (MFs) and 275 KEGG pathways (TableS S4 and S5). These DEGs were mainly involved in extracellular matrix (ECM) and immunity-related GO terms, such as BPs of ECM organization, T cell activation, and myeloid leukocyte migration, CCs of ECM, secretory granule membrane, and MFs of ECM structural constituent, immune receptor activity, cytokine binding, and cytokine receptor activity (Fig. 5B). KEGG analysis indicated the enrichment of immune-related pathways, such as cytokine-cytokine receptor interaction, cell adhesion molecules, focal adhesion, chemokine signaling pathway, and complement and coagulation cascades (Fig. 5C). Moreover, GSEA also found that ECM and immunity-related BPs and KEGG pathways were According to the median value of risk score, BLCA patients were divided into low-and high-risk groups (Fig. 7A). We found that patients in low-risk group had better survival than those in the high-risk group (Fig. 7B).
To evaluate the performance of the risk model, ROC curves were plotted and the AUCs for 1-, 3-, 5-year were 0.70, 0.72 and 0.71, respectively (Fig. 7C), indicating the high predictive accuracy of the risk model. Moreover, the risk model was validated in an external GSE13507 dataset (n = 165), and similar results were obtained (Fig. 7D-F), suggesting the stability and reliability of this risk model. In addition, we found that the risk score, age and N stage were independent prognostic factors in BLCA by multivariate regression analysis (Fig. 8A). To make it easier for clinicians to predict the survival probability of BLCA patients, a nomogram model was constructed based on risk score, age and N stage (Fig. 8B). The calibration curve was plotted to evaluate the nomogram's performance. Our data showed that the predicted probability was close to the actual survival (Fig. 8C)   www.nature.com/scientificreports/ The pyroptosis-related risk score was associated with immune-and chemotherapeutic efficacy. Previous studies have shown that the molecular subtype of BLCA can predict its clinical response to immune checkpoint blockade therapy and the neoadjuvant chemotherapy 2,4 . To further explore the potential capability of pyroptosis-related risk score, we firstly investigated the relationship between the risk score and immunotherapy efficacy. Our results showed that patients in low-risk group had lower TIDE score (Fig. 9A), indicating the better immunotherapy response for low-risk patients. As for chemotherapy, we focused on the small molecule compound, cisplatin, which was used as the first-line chemotherapy for BLCA. We found that the IC 50 of cisplatin were much higher in low-risk group (Fig. 9B), suggesting that patients in high-risk group were more sensitive to cisplatin treatment. Collectively, these results showed that our risk signature was associated with therapeutic efficacy of BLCA patients, and immunotherapy and neoadjuvant chemotherapy may be used, either alone or in combination, in the treatment of BLCA patients with high pyroptosis-related risk score.  www.nature.com/scientificreports/

Discussion
Due to the high recurrence rate, metastasis rate, and limited therapeutic agents, advanced BLCA remains a serious tumor disease of the urinary system in clinical practice 1,5 . Therefore, it is meaningful to develop novel diagnostic biomarkers, therapeutic targets, and risk model to predict patients' prognosis and guide chemotherapy and immunotherapy for BLCA. Numerous studies have highlighted the critical role of pyroptosis in modulating innate immunity and ani-tumor effects 13,20,21 . However, the specific molecular mechanism of pyroptosis remains to be fully elucidated in BLCA. Here, we integrated transcriptome data from TCGA and GEO database to comprehensively characterize pyroptosis-related clinical features and gene signature in BLCA. We identified two distinct subtypes based on two DEPRGs, namely, GSDMA and CHMP4C. Compared with patients in cluster 2, patients in cluster 1 had more serious clinicopathological features and worse OS. These two subtypes also displayed distinct characteristics of the TME. Furthermore, we identified 953 pyroptosis-related DEGs between two different subtypes, which significantly enriched ECM and immunity-related pathways. Based on these DEGs of the pyroptosis subtypes, we established a risk model to predict the prognosis for BLCA patients. By integrating the pyroptosis-related risk score and clinical variables, we constructed a nomogram to predict the overall survival of BLCA patients. Finally, we found that high risk score was significantly associated with high TIDE score and low IC 50 of cisplatin, demonstrating the clinical value of the prognostic model in predicting immunotherapeutic and chemotherapeutic efficacy. Our study not only provides clinicians with a practical prognostic model for stratification of patients, but will also further enhance our understanding of the molecular mechanism of pyroptosis in BLCA. Interestingly, while our manuscript was under review, several studies using similar strategy to develop PRGsbased risk models were published [36][37][38] . We performed a comprehensive analysis of every relevant aspect of the Gene  DOK7  IL9R  OVGP1  SPOCD1  HSD17B2  CTSE  FBP1  ISLR  CTHRC1  FN1  EFEMP1  COL6A1  DPYSL3  COL6A3  COL5A1  COL5A2  MFAP5  XAGE2  CCDC80  COL14A1  SCEL  TPM1  PDGFRB  ADAM12  FBN1  CYP26B1  NRP2  LOX  DPYSL2  PCOLCE2  MAP1B  ATP8B2  KCNE4  ADAMTS12  F2RL2  GAS7  MAP1A  ENPP1  KANK4  DSC1  STXBP6  ADAMTS16  PRKG1  GLI2  ALDH1L2  SLC19A3   www.nature.com/scientificreports/ risk model, including the pyroptosis-related gene sets, training and validation datasets, model construction methods and performance (Table S6). Overall, our risk model had better performance in predicting BLCA patients' prognosis. Moreover, these risk models had some slight differences in the number of PRGs, the strategy of model construction, and the training and validation datasets used. For example, Zhang et al. 38 used the gene expression patterns of the entire PRG gene set to conduct consensus clustering, while Chen et al. 36 , Deng et al. 37 and our study used DEPRGs. In our study, we initially used all PRGs to classify BLCA patients. All patients could be divided into two clusters ( Figure S2A, B). However, patients could not be clearly distinguished between the two clusters in the tSNE plot ( Figure S2C), and the overall survival of the two subtypes showed no significant difference (P = 0.44, Figure S2D). Therefore, we adopted the survival-associated DEPRGs to conduct the subtype analyses.
Our study identified significantly upregulated pyroptosis-related genes GSDMA and CHMP4C as favorable prognostic factors in BLCA patients. To date, the molecular functions of these two DEPRGs have not been well characterized in BLCA. GSDMA is one of the gasdermin family with pore-forming activity that functions as the key executioners of pyroptosis, which is mostly expressed in the epithelial cells of the skin, esophagus, and bladder and is silenced in the gastric cancer cell lines and tissues 11 . These evidences implied that GSDMA might be a tumor suppressor gene in gastric cancers. Saeki et al. also showed that restoration of GSDMA in human gastric cancer cell lines could increase their sensitivity to cellular apoptosis 39    Linear predictor www.nature.com/scientificreports/ that GSDMA was not expressed in normal colorectal tissues but was gradually overexpressed in matched cancer tissues, suggesting that it might act as tumor promotor gene in colorectal cancer 40 . A recent study showed that exotoxin B produced by Streptococcus pyogenes could cleave GSDMA and trigger keratinocyte pyroptosis, which verified the speculation that it acted as the host-defense protein in skin 41,42 . CHMP4C is an important member of the chromatin-modifying protein/charged multivesicular body protein (CHMP) family. These proteins are parts of endosomal sorting complex required for transport III (ESCRT-III), which is involved in degradation of surface receptor proteins and formation of endocytic multivesicular bodies (MVBs) 43,44 . During pyroptosis, ESCRTmediated cell membrane repair can negatively regulate pyroptosis triggered by the activation of GSDMD 45 . Chen et al. found that silencing of CHMP4C could obviously inhibit the proliferation, migration, and invasion of pancreatic adenocarcinoma cells 46 . Lin et al. also discovered that CHMP4C accelerated the cell proliferation and migration via activation of epithelial-mesenchymal transition pathway in cervical cancer 47 . These findings were not consistent with our data, in which CHMP4C acted as a tumor suppressor gene. More studies are needed to elucidate the reasons for the differences in the roles of CHMP4C across different cancers in the future.
In the past decade, breakthroughs have been made in the treatment of BLCA with immunotherapy. However, the heterogeneity of anti-tumor immunity is strongly linked cancer progression and response to immunotherapy 48 . TME is composed of tumor infiltrating immune cells (TIICs), blood vessels, fibroblasts, and ECM. Emerging evidence has also shown the significant roles of pyroptosis in modulating the TME, thereby affecting tumor development, progression, and therapeutic resistance 49 . In the present study, the pyroptosis subtype characterized by immune activation (cluster 2) was associated with a better prognosis. Immune checkpoint genes associated with immunosuppression, such as CD274, CTLA4, IDO1, and LGALS9, were markedly lower in patients in cluster 2. In addition, patients in cluster 2 had lower stromal/immune/ESTIMATE scores, which was consistent with the better prognosis. Accordingly, we also discovered that the relative abundance of 22 TIICs was significantly different between two pyroptosis-related subtypes, implying the key role of pyroptosis in shaping TME. Interestingly, we found that myeloid cells, including macrophage M0/M1/M2 and activated DCs, showed significant differences in the two subtypes. Cluster 2 with better prognosis showed higher infiltration of activated DCs, suggesting that they play a positive role in BLCA. Although the fraction of DCs in the TME was small, they had the potent capability to foster T cell immunity and enhance immunotherapy responses 50 . Tumor-associated macrophages are polarized into two main phenotypes, namely M1 macrophages and M2 macrophages, which inhibit or promote cancer progression, respectively. In this study, the fraction of M2 macrophages with immunosuppressive phenotypes in TME was significantly lower in cluster 2, which was associated with the favorable prognosis. Intriguingly, we also found decreased infiltration of M1 macrophages in cluster 2. We speculated that this might be related to the stronger plasticity of macrophages [51][52][53] . To improve the accuracy of dissecting the TME composition, sing-cell transcriptome datasets might be required to analyze the proportion of macrophage under different polarized states in the future.
We further explored the relevant signaling pathways and functions of DEGs between different pyroptosisrelated subtypes. Both KEGG and GO enrichment analysis indicated that DEGs significantly enrich a variety of biological functions related to ECM and immunity, including ECM organization, T cell activation, myeloid leukocyte migration, and cytokine-cytokine receptor interaction. These findings were consistent with the estimated stromal and immune score by ESTIMATE algorithm. In addition, we have also established a nomogram for predicting prognosis at 1, 3, and 5 years, which can help clinicians to accurate stratification of BLCA patients. Last but not least, we found that this prognostic model could guide immunotherapy and chemotherapy, which would bring clinical benefits to BLCA patients.
There were some limitations for this study. Firstly, our findings were solely based on public datasets that were composed of retrospective data. Large-scale prospective clinical cohort validation in the real world is needed. Secondly, we used a relatively small number of paracancerous tissues (n = 19) in the TCGA-BLCA dataset as normal controls to identify DEGs compared to the BLCA tissues (n = 414), which may cause bias. In addition, paracancerous tissues, the tissue adjacent to a solid tumor, although histologically normal, are not normal tissues after all, which might suffer from a field effect. Thirdly, the molecular mechanism of pyroptosis shaping TME characteristics in BLCA was still unknown. Additional in vivo and in vitro functional studies should be investigated experimentally in the future.
In summary, we utilized a comprehensive bioinformatics approach to dissect the molecular mechanism of pyroptosis in affecting TME composition, clinicopathological features, and patients' prognosis in BLCA. We also constructed a pyroptosis-related prognostic model to inform immunotherapy and chemotherapy. Our findings not only provided novel insights into the role of pyroptosis in BLCA, but also supplied potential independent prognostic factors for diagnosis and therapy of patients with BLCA.